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For the exactly solved reduced BCS model an electrostatic analogy exists; in particular it served 
to obtain the exact thermodynamic limit of the model from the Richardson Bethe ansatz equations. 
We present an electrostatic analogy for a wider class of integrable Hamiltonians with pairing force 
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the analytical results, we compare them with numerical solutions of the Bethe ansatz equations for 
finite systems at half-filling for the ground state. 
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INTRODUCTION 



Pairing force interactions have been successfully em- 
ployed to explain phenomena in differenL contexts such 
as superconductivityEI, nuclear physicsS, QCDd and 
astrophysicso. The original idea traces back to Bardeen, 
Cooper, and SchriefFer (BCS) who proposed pairing of 
time-reversed eleetrons as the crucial mechanism for 
superconductivityEI. Physical implications of the corre- 
sponding Hamiltonian have been extracted resorting to a 
variety of analytical or numerical techniques. The mean 
field ground state was shown to be exact in the limit 
of large uumher of electrons where fluctuations can be 
neglectecu'ElQ. This constituted a great success of the 
BCS variational ansatz. However there are relevant phys- 
ical situations where approximations are not reliable and 
exact treatments are highly desirable. Fortunately, by 
properly choosing the pairing couplings, the model ad- 
mits an exact solutions. 

A simplified, but still non-trivial model is the "reduced 
BCS model" (BCS model, for brevity) which assumes 
a uniform pairing g among all the electrons within the 
Dcbye shell. This model is integrablelll and was diago- 
nalized long ago by Richardson through Bethe Ansatz 
(BA)la (see also Ref. Only very recently this gener- 
ated a lot of interest both in nuclear and condensed mat- 
ter physics whose communities benefited from the simple 
algorithm of Richardson's BA soluiip*i|-t«-4ackle the BCS 
model in the canonical ensemblcli2lllirli30. Much work 
has been done to merge the model in the schemes of the 
Quantum Inverse Scattering (QIS) and Conformal Field 
Theory (CFT), which are modern arenas where quantum 
integrability and exact solutions can be treated on equal 
footing. These studies allowed significant steps forward. 
QIS studies identified the Richardson BA solution as the 
quasi classical limit of the exact solution of (tjwisted) dis- 
ordered six vertex models of the XXX-typeEJIlj. This 
paved the way towards the exact |^valuation of corre- 
lation functions of the BCS modelEJO. The field the- 
oretical study of the BCS model is due to Sierra and 
coworkers. The first step was to relatcr-tjie, Richardson 
BA solution with WZNW-su(2)fe modelfJlal3. This stim- 



ulated the discovery that the BA solution is the quasi 
classical limit of the Babujian's off ^b^X BA of the (un- 
twisted) disordered six-vertex modelEJo. The field the- 
oretical origin of the integrability of the BCS model has 
been clarified definitelytiJ to be a twisted Chern-Simons 
(CS) theory on a torus. The Richardson BA solution 
together with the underlying integrability of the theory 
arises from the Knizhnik-Zamolodchikov-Bernard equa- 
tions (which are the Knizhnik-Zamolodchikov equations 
on the torus). The emergence of the CS theory is quite 
interesting and relates at a formal level the BCS the- 
ory of superconductivity with the Fractional Quantum 
Hall Effect (FQHE)ll3 (also pointing out important dif- 
ferences). In fact, the exact BCS wave function admits 
a "Coulomb gas" representation, which corresponds to 
the "Coulomb gas" representation of the Laughlin wave 
functions (plasma analogy )lla'E3. Accordingly, also for 
the BCS model an electrostatio. analogy does exist. It 
was illustrated first by GaudinE3: the Cooper pair ener- 
gies are obtained as equilibrium positions of N mobile 
charges in the background of Vl fixed charges and a uni- 
form electric field of strength iJg. This analogy was used 
by RichardsonO and GaudinCJ to obtain the thermody- 
namic limit of the Richardson BA equations in two dif- 
ferent approaches. Both reproduced the BCS mean-field 
gap enuation, hence confirming the statement of Bogoli- 
ubovO that at T = the mean-field results become exact. 
Recently, Sierra and coworkers revived the attention on 
the thermodynamic limit of the BCS model, comparing 
Gaudin's long forgotten results with the numerical solu- 
tion of the pRichardson BA equations for a finite number 
of particlesEj. 

By generalizing the integrals of motion of the BCS 
model, the class of known integrable pairing-force Hamil- 
tonians could be enlarged corisiderably towards models 
with non-uniform interactionsESEZI (for bosonic versions 
see the Ref.^). Here, the term "non-uniform" indicates 
that the interactions depend on the energy levels occu- 
pied by the interacting electrons. As the uniform BCS 
model, also these models emerge from the quasi-classical 
limit of the QIS approach for, twisted disordered six ver- 
tex models of thepfXZ-typdlj'Ej. The untwisted case, 
studied in the Ref.EJ, corresponds to large effective pair- 



2 



ing interaction. 

In the present work we generaUze the 2c?-electrostatic 
analogy to the class of integrable Hamiltonians obtained 
in Ref. 26 , The thermodynamic limit of these Hamiltoni- 



ans is obtained. We compare the obtained analytical re- 
sults with numerical solutions of the BA equations. The 
electrostatic analogy of the Richardson-Sherman equa- 
tions is a screening condition for a total electric field pro- 
duced by charges distributed in two-dimensional space. 
Accordingly, the electric field obeys a Riccati-type differ- 
ential equations. 

The paper is laid out as follows. In the next section 
we review the integrability of pairing force Hamiltoni- 
ans that arises from the underlying infinite dimensional 
Gaudin algebra ^/[s^(2)]. In section HI we present the 
electrostatic analogy, review basic facts of 2(i-dimensional 
electrostatics and apply them to obtain the thermody- 
namic limit of integrable Hamiltonians with non-uniform 
pairing couplings. At the end of the section, wc obtain 
Gaudin and Richardson's results for the BCS model as a 
limiting case of our equations. Section IV is devoted to 
conclusions. In Appendix A we sketch some mathemati- 
cal aspects connected with the integrability of the mod- 
els. In Appendix B we present the connection with the 
Riccati equations by reviewing the work of RichardsorO. 
In Appendix C we collect details of the calculations. In 
Appendix D we discuss some features of the ground state 
and possible routes towards the study of the excitations. 



II. EXACTLY SOLVABLE PAIRING MODELS 

The Hamiltonian for N charged particles interacting 
through pairing force 



has been provedcS to be exactly solvable if the Coulomb 
Uij and the pairing strength gij are restricted to have the 
following form 



9ij 



pg— 



for i =/= j 



smh[p(ui — Uj)] 

Wij = A - pg{£i - £j) coth[p(ui - Uj)] 

(2) 

where|-et are the single-particle levels, while for conve- 
nienceLJ we fix gu = (j), Uu = A — (j), where 20 = 
pgJ2i^ji^i ~ ^j) coth [p{ui — Uj)]. The quantities A, u, 
and Ui are arbitrary real parameters while p can be real or 
pure imaginary. To study the integrability of the Hamil- 
tonian (|l|) it is convenient to reformulate it as a spin chain 
model. In fact it can be expressed (up to a constant term) 
as 

- \ E^,=i <?y- {K+KJ + KrK+), (3) 



where K+ = c]^c]^, Kj - ,K-^ = 

\{nj-\ -\- Uji — 1) , are su{2) operators. The Hamiltonian 
(|l|) acts non trivially only on the Hilbert space Ti of dou- 
bly occupied and empty level pairs; the singly occupied 
levels are excluded from the dynamics associated to (|^). 
This phenomenon is commonly called the "blocking of 
singly occupied levels" . It reduces the problem of find- 
ing the spectrum of H in the restricted Hilbert space 
TL and it is an important ingredient for the integrabil- 
ity of the Hamiltonian (|l|). Ultimately the latter resides 
in the algebraic connection of the pairing Hamiltonians 
with the infinite dimensional Gaudin algebra (see Ap- 
pendix A) g[sl{2)] = span{K^{u),K''{u) : w G C}: the 
diagonalization of the Hamiltonian is equivalent to the 
diagonalization of the correspondent of Casimir operator 
in 5[s^(2)], K{u) . This is achieved either through al- 
gebraic or coordinate-wise BA, and corresponds to the 
quasi-classical limit of the inhomogeneous XXZ model. 
The set of operators in involution are the residues oiK{u) 
and the Hamiltonian is a polynomial of these operators 
which are, by construction, the integrals of motion of the 
model. 

In the following we will see these statements at work for 
the model (|), (|). Then we will recover the BCS model 
as a limit of the pairing couplings. 
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FIG. 1: Plot of the pairing and Coulomb couplings in Eqs.(^) 
with Ei = i, Ui = Ei, A — 2.5 and g — 1. In the numerics done 
in this work these couplings have been considered for various 
values of the rescaled g. 



1. Non uniform couplings. 

We review the integrability of the pairing models (|^). 
The key observation is that the integrals of motion fj of 
the model are residues of t(u), defined in Appendix A, 
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in u = Ui. In fact, H can be written as 
n n 



H - ^2e,f,+A^f,f,-^0,K2, (4) 



h = iff + , 



2 sinh — Uj)] 



(5) 
(6) 



+K-Kj coth - Mj)] 



The 5 — > oo-limit of this Hamiltonian was found in 
Ref. |30[ From Eq.(|) it is evident that [H,fj] = 0. By 
these integrals of motion, the model is connected with, 
the anisotropic Gaudin Hamiltonians 5; . The propertjO 
f;] = for j, Z = 1, . . . , comes from [t(u), t(u)] = 
E^. The common eigenstates of H and ti are 



becomes connected with isotropic Gaudin Hamiltonians 
Sj. In fact, H can be written as 

o n 
Hbcs = ^ + g ^ Tj + const. (12) 

i— 1 — 1 

= Kf + Ei , 



-9 



(13) 
(14) 



where the spin vectors are: Kj := {K^ , Kj , K^); = 
± iKj. The eigenstates of both the BCS Hamilto- 



nianHH and its constants of motion Tj-Q are 



N 



(15) 



AT 



1*) = l[K+{e„)\0) 



(7) The eigenvalues of the Hamiltonian and of the integrals 
of the motion are respectively 



where |0) is the Fock vacuum and N is the number of 
pairs. The eigenvalues of the Hamiltonian H and of the 
constants of the motion are 



N 



n N 



£ = ^pg ^ Ej coth [p{ea ~ Uj)] + AN'^ (8) 



£bcs — 2 ( 



^ 1 



(16) 



1 1 
E:r^-i (17) 



n = -| y^pcoth[p(ea - Ui)] + |y^pcoth[p(Mj- -m)]--. 
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2 Cq, Ei 4 ^.^^ Ej Ei 

where Ca are solutions of the Richardson- Sherman (RS) 
equations 



(9) 

The quantities Cq are solutions of the following set of 
equations 

2 ^ a 

- + 2^pcoth [p{ef3 — Ca)] — ^pcoth [p{uj — Cq)] = . 



I 1/2 ^ 1 

- + > > 

9 ea - £j e„ - ep 



(18) 



/3 = 1 



which are the BA equations of the model (pi) 



(10) 



(the parameters Cq, are one half af the parameters Ea 
originally defined by RichardsonO) . The eigenstates, 
eigenvalues and the RS equations can be obtained in 
the p ^ limit of the BA equations ( [l0| ) with the 
identification Ui — Ei, g = g. By the same limit the 
anisotropic (XXZ) Gaudin model reduces to its isotropic 
version (XXX). 



8. Uniform couplings: the BCS model 

We now review the exact solution together with the in- 
tegrability of the pairing Hamiltonian with uniform cou- 
pling constants. In this case the Hamiltonian(|l|) reads 

o n 
Hbcs = ^ ^^^^^ " f E i^^^^i + ^t) 

i—l — l 

H can be dkectly diagonalized through coordinate-wisei 
or algebraicCil BA. Nevertheless its diagonalization can be 
achieved along the lines depicted above. The integrals 
of motion of the BCS model are the residue of t{u) 
in w = 2ei. By these integrals of motion, the model 



III. ELECTROSTATIC ANALOGY FOR 
EXACTLY SOLVABLE PAIRING MODELS 

The RS equations ( p^ ) are stating that the total force 
acting on the, unit charges located in Cc, due to the 2d 
electric fieldO generated by the Vl charges q/ixed = —1/2 
fixed in Si, by the remaining — 1 mobile unit charges 
Qmobiie = +1 and by the external constant electric field 
of strength ~l/g is zero. In other words, the solutions 
of the RS equations correspond to the (unstable) equi- 
librium configurations of N charges in the complex plane 
under the influence of a given electric field. T|ljiis electro- 
static analogy was first pointed out by Gaudinu. It allows 
the exact access to the thermodynamic limit of the BCS 
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modelii by exploiting the formalism of complex analy- 
sis. The analytic structure of the ansatz electric field is 
prescribed by the positions of the charges. 

We will now present an electrostatic analogy for the 
generahzed BCS models (|l|), (|)- 

In the variables 



tanh (pea) 



the Eqs. (no) are algebraic 



tanh(p7ii) 



(0) 

9+ 



(0) 



l/p Za + l/p 



n 



1/2 



N 

E 



(19) 



0, 



(20) 

where 2g£°^ ^ ^l/pg-[n-2{N -l)]/2. We shall assume 
that the parameters ui are given by Ui = f{ei) where / 
is a monotonic function. For the case p = i note that 
due to the 7r-periodicity of the transformation (|l^) the 
Ui can be restricted to lie in (— 7r/2, 7r/2). Hence, the 
Bethe equations for the generalized BCS models (10) can 
be recast in the same form as the original RS equations 
Eq. (|l^), except that the constant external electric field 
of strength — l/g in the homogeneous case is replaced by 

the pr|eaence of two isolated charges —q± fixed at points 
±l/p.E£l Note that for p = i the charges are complex. 
The interpretation of this is discussed later on. 



A. Basics of 2d electrostatics 

We sketch briefly the main ingredients of two- 
dimensional electrostatics in terms of holomorphic func- 
tions in C (i.e. real harmonic functions in IR'^). 

Let the electric field E = [E^, Ey) be associated to the 
complex number E = E^ — i Ey. The Maxwell equations 
can then be summarized as 



d,E = ^{divE + i (Curl^),) = 7r(p + i ijmag)z) (21) 

with real charge density p and real magnetic monopole 
current density jmag in z-direction; the derivatives are 
dz = {dx — idy)/2, dz — {dx + idy)/2. Any integral over 
a closed curve T gives 



1 

2711 



[dzE{z) = — / {Exdy ~ Eydx) ~ i{Exdx + Eydy) 
Jr 27r Jr 



= 1. lE-ds-1- 



2n 



2tt 



E dl 



i/r 



where Qr and /p are the total electric charge and mag- 
netic monopole current enclosed by T (in the following 
we refer to Qr — i/r as charge) . The surface element ds 
is a vector perpendicularly pointing outwards of T with 
the same length as the line element dl. 

A charge q contributes to the electric field with a simple 
pole with residue q. A line of charges gives a holomorphic 



function on a Riemann surface with a branch cut along 
the charge line. The discontinuity of the field in cross- 
ing the cut gives the charge density ^ {E^{z) — i?+(z)), 
where E-{z) {E+{z)) are the limiting values of the field 
when z tends to the cut from the right (left) with respect 
to the orientation of the curve. 



B. Thermodynamic limit 

The thermodynamic limit of the model (|^) can be ob- 
tained in the following way: we first divide Eq. (|2^) by 
and define the (positive) charge densities 



a{za) := 



1/2 



fl{Xj + l - Xj) 

1 



To obtain a sensible thermodynamic limit, we assume i) 
that the pairing strength scales as g = G/fl, with fixed 
G; ii) that the Debye shell defined by the end points ei, 
en does not depend on fl; iii) that the number of pairs 
increases with fl according to = vfl, where v is the 
filling. In the limit — > oo. Equations (|2y) then become 



i(o) 



l/p 



l/p 



dx 



p{x) 



dz 



(22) 

where the integrals are meant in the sense of the prin- 
cipal value V. After the transformations ( |l9| ) the trans- 
formed Debye shell L is still a segment of the real axis 
having end points qq = tanh {p ui)/p, ba = tanh (p uq)/p\ 

the isolated charges are — Q^"* — — 1/2 ± l/(pG)); 

the density p is determined by the single-particle energy 
density as p[x) = p^{e{x))/[f'{e{x))[l —p^x^)], where 
e{x) ~ /^^(Arctanh(pa;)/p), and it fulfills 



deps{e) = dx p{x) = 1/2 



The curve F and the density a have to be determined, 
with the constraint 



\dz\a{z) 



The idea of Gaudin was to construct the total elec- 
tric field E{z) as a function with an analytic structure 
prescribed by the actual charge distribution. He fur- 
ther assumed that eventual solutions of the equation 
( p^ ) are arranged in K piece- wise differentiable arcs: 
F = Fp U [Jn=i^n, where Fp is the common support 
of p{x) and cr(x) which is not contained in the arcs 

Tare = [Jn=i^ri- For example, a charge distribution 
along the line [a, b] in the complex plane with total charge 
q and constant charge density would lead to the total 
electric field E{z) = {q/\a — b\) ln[(z — b)/{z — a)] which 
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diverges at the end points of [a, b]. Since in the case un- 
der consideration the electric field has to vanish at the 
end points of F, we need the charge density to vanish 
"sufficiently fast" there. An example of such a function 
could be any [{z — a){z~ b)]°' In with any a £ ^IN or 

just [{z — a){z — b)]~ with odd a. Both functions indeed 
vanish in a and b and have a branch cut along [a,b]. The 
choice of the admissible functions is finally restricted by 
imposingpthat the free charges are distributed along a fi- 
nite set rE3'Ej. In fact, for a ^ i]N we either have infinite 
branch cuts or divergences at the end points. We want 
to stress that the run of the branch cut can be chosen 
in an arbitrary way as long as it joins continuously the 
branch points a and b. The run is fixed requiring that 
it is an equipotential curve of E{z). The ansatz function 
for the total electric field is 



we have that 



E{z)=S{z) 

K 



Q- 



z — \/p z + l/p 
S{z) = \/{z - an){z - bn). 



dx 



(23) 



n=l 



Note that the sign of the complex square root is uniquely 
given, once Tare is fixed, by yj{z ~ an){z - bn) = 
yj\z - an\\z - 6„| exp ^ argp^(z - a„)(z - 6„). The ar- 
gument function argp {z — an){z — bn) is the sum of the 
angles between z G C and the two end points of r„ , which 
has a 27r-discontinuity crossing r„. The unknown quan- 
tities are the end points (a„, 6„) of Tare, Q± and (p{x). 
To give a physical interpretation of such an ansatz func- 
tion, look at S{z)Q± and S{z)(f{x) as screened charges 
and charge density, respectively. From the symmetry of 
the distribution of the fixed charges on L and in ±l/p 
it follows that if a„ has non-zero imaginary part, then 
bn = o,n, and we choose the orientation of r„ from a„ to 
bn, with Im(a„) < 0. If a„ is real, so must be bn- We ar- 
gue however that this occurs only for distinct "critical" 
values of the pairing coupling at which the imaginary 
part of complex conjugate a„ and 6„ vanishes and hence 
o-n = bn. The corresponding arc is then describing a 
closed curve. Therefore we consider S{z) to don't have 
cuts along the transformed Debye shell L (see Appendix 
We emphasize that this does not mean that a{z) has 
no support on L. A non-zero solution charge density on 
L is accounted for by the screening density S{x)(p{x). 

After having determined the a„,6„, the curves r„ are 
found as equipotential curves of the total electric field 
Eq. (p3|) . The density a is then determined by the dis- 
continuity of E in crossing T. 

In the following, all the unknown quantities are ob- 
tained following a procedure d la Gaudin on which we 
report in the Appendix We will calculate the complex 
integrals involved in this procedure by exploiting basic 
knowledge of electrostatics. 

Since the screened charges must tend to their bare val- 
ues when the point z is sufficiently close to the source. 



S{±l/p) Qi 



(24a) 



In a similar way, using (j{z) — p{z) = [E^{z) 
E+{z)) / {2T:i) , we find 



S{x) (p{x) — cr{x) — p{x). 



(24b) 



Next, imposing that the electric field asymptotically 
goes as E{z) ~ l/z^ (this is the leading term in the 
multi-pole expansion, since the total charge is zero), we 
obtain the following equations 



dxx'^'ifix) = (Q+ + (-l)"Q_), 0<n<i^.(25) 

L P" 



We present a detailed derivation of these equations in 
Appendix There we also show explicitly that the 
condition ^^a{z)\dz\ = v is automatically fulfilled with 
Eqs.dl). 

Using Eqs.(^) we can express Q± and f{x) in terms 
of S{±.l/p) and S{x) respectively. The r„ are finally 
obta^^as equipotential curves of the total electric field 



Rc dz' E^ {z') = J E.^dx + Eydy = , 



(26) 



The density cr is determined by the discontinuity of E 
along r, according to: 



<j{z)^p{z) + -\E^{z) 

TT 



(27) 



It is worth noting that Eqs.(p5|) are a set of if + 1 real 
equations for the 2K real parameters determining the 
end points of r„. Gaudin stated that seemingly in any 
case there is a finite number of solutions. From a phys- 
ical perspective, though, we expect that the K — 1 free 
parameters span a family of curves corresponding to a 
band of excitations of the systeHi. This conjecture will 
be verified in forthcoming worl£3. 

We summarize all the conditions found for the un- 
knowns a„, bn, ^, Q 
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K 



S{z) = Y]_ Vi^ - 0.n){z - bn) 



E{z) = S{z) 



Q- 



z — 1/p z + 1/p 

v-l/2±l/pG 
2S{±l/p) 

a{x) - p{x) 
Six) 

a{z) = p{z) + -\E^{z)\ 

TT 

The arcs are then determined by 

,,p{x)~cj{x) _ Q+ + {-rQ- 



dx 



ip{x) 



Q± 

fix) 



dxx^' 



Six) 



p" 



n^O,...,K 



Z 

r„ : Re Jdz' E- (z') = . 



As we stated above, the density p of the transformed 
variables x — tanh {pf(e))/p is connected with the given 
density p^ of the single-particle levels via 

p(x)=Pe(e(x))/[/'(£(a;))(l-pV)], 
e(x) ^ /^"^(Arctanh(pa;)/p) . 

The energy of the state is finally given by 



£ — av^ — G ldxp{x)e{x) \dz\a{z) 



1 — p xz 



av^+p^Gv i de p£{e) e x{e)+ 



f'i^ix)) Jr z-x 



(28) 



(29) 



where a = AVL. 



1. Comparison with the thermodynamic limit of the BCS 
model 



In the limit p ^ 0, the two isolated charges Q^'' take 
an infinite value and they are displaced to infinity in such 
a way to give a uniform electric field 1/G. In this limit, 



the charges Q± behave as 



K 



11/ 11 

2G^ 2 l^" 2^2G 



K \ - 

^(On +bn)\p 
n=l ) . 



K-1 



where we kept into account the dependence of the rela- 
tive determination of S{±l/p) upon the number of arcs. 
Thus, the ansatz field reduces to 



Ep^oiz) = ^o(^) = S{z) dx 



(p{x) 



(30) 



Eq.(p4b|) still hold, while Eqs.(|25D simplify to: 



dxx''ip{x) =0, 0<n<K-2, (31a) 

/ dxx^-V(^) = (31b) 

f 11^ 

/ dxx'^ifiix) ^ — ^{an + bn), (31c) 



The last two equations come from the fact that the field 
should go asymptotically as Eo{z) ^ — l/G + Qt/z in 
this limit, since the total charge is now Qt = z/ — 1/2. In 
the approach by Gaudin the last condition was obtained 
by imposing J\dz\a{z) — v. This extracts the residue 
at infinity, which has to be such that the total charge is 
zero. In our approach, we already started from a globally 
neutral system and hence this normalization condition for 
a{z) was automatically fulfilled (see Appendix |c|) 

We compare our findings for integrable pairing models 
with the results for the BCS model obtained in Refs. 0, 
p3 and reconsidered recently by Sierra and coworkersEj 
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K 



coincide if the pairing coupling strength u is zero) 



S{z) = ]^ v^(z - a„)(z - bn) , 

71=1 

ip{x) 



Eo{z) = S{z) dx 



Z — X 

, , a{x) — p{x) . . , . 1 I ^ , X 

'^^■'^^ " S{x) ' '^(^) + 

The arcs are then determined by 



I 



dx : 



dx 



S{x) 



= , n = 0,...,K -2 



„K 



S{x) G 



K 



J Six) - 2G + - - - ^ 



Z 

r„ : Re Jdz' Eo-iz') = . 




G = 0.50 
G= 1.00 
G= 1.50 
G= 1.75 
G= 1.95 
G = 2.00 
charge 
L 



FIG. 2: Plot of the pairing parameters z in the complex plane 
for the hyperbolic model corresponding to p = 1 in Eq.(|^. 
The symbols refer to the numerical solutions of Eqs.(^) for 
fl = 200, A'^ — 100, Ui = Ei = 0, . . . ,2 and uniform level spac- 
ing. Different colors and symbols represent different values of 
G, as shown in the legend box. The lines are the plots of the 
curves determined by Eq.(p^). If we choose the Debye shell 
to be (—1, 1) as in the trigonometric case, the arc for G — 2 
closes at infinity. 



2. Numerics and discussion of the solutions 

To solve the BA equations (J^, @ we have chosen 
the parameters Ui = £i at half filling = 2N . The corre- 
sponding interactions Uij and gij are shovifn in Fig.|l]. In 
the trigonometric case the Debye shell ranges in (—1,1); 
in the hyperbolic case it ranges in (0, 2) for numerical 
convenience (see the caption of Fig.^). The ground state 
of the system (which is the only state we consider in the 
numerics) is obtained by evolving the G = state 



lim Za 



ae (l,...,iV). 



(32) 



to a finite value of G. At certain value of G the equations 
are singular since a couple of pairing parameters Z2\~i 
and Z2X coincide with the energy level X2\-i. Such a sin- 
gular behaviour can be smoothed out using the procedure 
developed by Richardson in Ref. ^ In the following we 
briefly summarize it. The corresponding divergences are 
removed by the following transformations 



Z2\~i = A\ - iBx, Z2\^ A\ + iBx, (33) 



and then 



. X2\-l + X2\ 

Xx — Ax , Yx 



Bl 



(34) 



where Sx = {x2x ~X2X-i)/'2 is half the energy spacing be- 
tween the corresponding energy values (with whom they 



By the first transformation Eqs. ( pO| ) can be written 
in a form which is manifestly real, whereas the second 
transformation ( ^^ removes the divergences. 

Now we discuss some features of the arcs shown in 
Figs.|^ and For K = 1 (we recall that straight lines 
on the real axis are not counted by the index n) , the end 
points are uniquely determined (if theA!-|exist) and the 
curve F corresponds to the ground stata23 of the system. 
In this case the end points of the arc, A ± iA, are deter- 
mined by the two coupled equations (obtained as linear 
combinations of Eqs. 



dx 



dx 



(1- 


■\-px)p{x) 




- A)2 A2 


(1- 


-px)p(x) 




- A)2 + A2 



= 2Q+ = 
= 2Q_ = 



l/pG + iy-1/2 
V(A - l/p)2 + 
-1/pG + i^ -1/2 

V(A + i/py +^ 



Why this set of equation corresponds to the ground state 
is discussed in some detail in Appendix 

In the trigonometric case (p = i ) , the equations are 
both complex, and they are conjugate. Thus, in order 
to find A and A it suffices to solve separately real and 
imaginary part of one of them. In the hyperbolic case 
(p — 1), assuming a uniform single-particle energy den- 
sity p^{e{x)) = po (i.e. p{x) — po/(l — x^)) and a De- 
bye shell spanning the interval [0, 2ujd], we have that the 
system above admits a real solution (and hence two co- 



inciding end-points) 
A =0 

A = sinh(2tj£))/ (cosh2w_D — exp [(i/ — l/2)/po]) , 

if 1/ {poG) — 2uj]j . This corresponds to the external curve 
in Fig. (1^ 




Rez 

FIG. 3: Plot of the pairing parameters z in the complex 
plane for the trigonometric model which correspond to p = i 
in E q.(B[ ). The symbols refer to the numerical solutions of 
Eqs.(|^for = 200, N = ICQ, u, ^ Si = -1, . . . ,1 and uni- 
form level spacing. Different colors and symbols correspond 
to different values of G, as shown in the legend box. The lines 
are the plots of the curves determined by Eq.(^6|) . 

Situations corresponding to ii' > 1 lead to excited 
states. We point out that the independent equations 
( p5| ) are in number oi K + 1 leaving undetermined K ~ 1 
parameters. This might be an evidence of "bands" of 
excitations above the ground stata2a. 

Looking at the behaviour of the quasi-momenta solu- 
tions of the BA equations, we can note that the more 
intense is the pairing constant G the more evident is the 
tendency of the quasi-energies to be complex. In Fig. ^ 
we compare this tendency for the hyperbolic, trigonomet- 
ric, and uniform BCS models. We found that the pairing 
tendency increases from the trigonometric over the uni- 
form to the hyperbolic model. Looking at the explicit 
couplings plotted in Fig. 0, we should conclude that the 
pairing and Coulomb interaction in Eqs. are compet- 
itive couplings. 



IV. SUMMARY AND CONCLUSIONS 

We have discussed the electrostatic analogy for the 
exactly solvable pairing force Hamiltonians found in 
Ref. These models generalize the BCS model in that 
the coupling constants are not uniform (in the space of 
quantum labels see Eqs. (|l|)). The ordinary BCS 



e 



O Unifoim 
□ Hyperbolic 
O Trigonometric 



ooooooo^°o° 




-0.6 

Ree 



FIG. 4: For an exemplary value of G = 2 the solutions 
for p — l,i corresponding to hyperbolic and trigonometric 
models respectively are compared with the solutions of the 
BA equation for the uniform case p = which are the RS 
equations (the Debye shell ranges in ( — 1, 1) for all the three 
cases) . 



model is obtainable as a limitil p ^ of the models 
(|l|) , (|^) . The electrostatic analogy exists because the B A 
equations for the "rapidities" Cq can be recast in a spe- 
cial algebraic form. For the class of models we deal with 
such an algebraic form is obtained by the change of vari- 
ables: Za = tanh(peQ), Xi = tanh(pei). After this trans- 
formation the Bethe equations express the condition for 
an equilibrium alignment of mobile charges qmobUe = +1 
in the 2d-real or complex plane in a neutralizing back- 
ground of q fixed — —1/2 fixed charges placed in the po- 
sitions Xi and two further charges in the positions ±l/p. 
Remarkably, this analogy is very effective to obtain the 
exact thermodynamic limit (i.e. the large N limit) of 
the models. In this limit the mobile charges (which are 
the solutions of the Bethe equations) arrange themselves 
along equipotcntial curves F = IJ^^ F„ of the total elec- 
tric field with vanishing charge density at its extremi- 
ties. The electrostatic analogy presented in this work 
is by no means unique; the charge distribution can be 
transformed by the whole group of conformal Moebius 
transformations. We have obtained explicitly the ther- 
modynamic limit for the ground state configuration at 
half-filling (see Fig. (||), (||)). We note, however, that 
the approach is also valid for excited states of the sys- 
tem, though the explicit calculation can involve technical 
subtleties (discussed in Appendix ^ . 
The electrostatic analogy together with the construc- 
tive equations for the BCS model found by Gaudin and 
Richardson is demonstrated to be obtained from the limit 
p — !■ of the results presented here. A comparison be- 
tween the hyperbolic, trigonometric, and uniform cases 
is seen in Fig. (^. Since this represent an inversion of 
the behaviour one would have expected from Fig. (^ we 
conclude that the pairing and Coulomb interactions in 
Eqs. (0) are competitive couplings. The exact thermody- 
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namics of integrable pairing models is one of the future 
goals. 
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We note that RS equations ( |A3[ ) are intimately related 
to the algebraic structure of ^[sZ(2)] since they act as 
constraints on the lowest weight k{Ea)- The difference 
between the BCS and Gaudin model results in a different 
constraint imposed on the lowest weight vector of t/[sZ(2)] 
which leads to different sets £, £' of solutions of t he BA 
equations {£' is spanned by the solutions of ( |A3| ) when 
g — > oo). This fact has been used to extend the Skljianin 
theorem for the Gaudin models to the BCS modeO. 



APPENDIX A: THE GAUDIN ALGEBRA 



2. Hyperbolic/Trigonometric Q[sl{2)] 



The Gaudin algebra 5[sZ(2)] is constructed from sl{2). 
The sl{2) "lowest" weight module is generated by the vac- 
uum vector \0)j, K~\0)j = , Kj\0)j = kj\0)j where 
kj is the "lowest" weight {kj = — 1/2 for spin 1/2 which 
is the case considered for electrons). The infinite dimen- 
sional Q[sl{2)] is generated by 



In this case: :— p/sinh[p£_], -0(0 pcoth[p^] , 

k(^) EE fc(^) = Y^f^ipkj coth[p{^ - Uj)]and A = g. The 
operators ( |Ai| ) obey 



K^{v) — cosh[p(ti — w)]K^{w) 
sinh[p(w — w)] 



n 



(Al) 

where ^ e C and Uj e IR. The module of ^?[s/(2)] is char- 
acterized by the vacuum |0) = 'Sif^i\0)j: K~{^)\0) — 

, K^mO) = <m , where k{0 := E'Li 



iK-i.),K^.)^^2pEi^^, 
smh[p(v ~ w)\ 

where v ^ w e C. 

The BA equations for the corresponding Hamiltonian 
(IH) are 



is the lowest weight of g[sl{2)]. The element of fc(e^) ^ _ + ^ pcoth[p(e/3-e„)] , a = l,...,N . (A4) 



Q[sl{2)] corresponding to the su{2) Casimir operator is 

The generating function of the integral of the motion of 
the pairing models is related to K{(^) through 



TiC) = -2AK{0 + K^0+cm 



(A2) 



where c(^) is a C-function and A is a real parameter. The 
property [t{v),t{w)] = arises from the quasi-classical 
limit of the sl(2) QIS theory; this is the ultimate reason 
for the integrability of the BCS model. 

In the present paper (j) and ip are either rational or 
trigonometric /hyperbolic functions. 



Rational g[sl{2)] 



In this case: 0(0 = VXO VC, <0 = HO 



J2i=i kj/{£, — Uj), and A = g. The operators ( |Al| ) obey 



[K-{v),K+{w)] 



V — w 
K'-{v)-K%w) 



where v ^ w G C. The Bethe equations are the 
Richardson-Sherman (RS) equations 



1 ^ 1 

^0(60) = - + X! — I 



Also in this case the BA equations act as constraints on 
the lowest weight k[ea) of the trigonometric t/[sZ(2)]. 



APPENDIX B: THE RICHARDSON ROUTE AND 
THE RICCATI EQUATION 

In this appendix we first emphasize the key role played 
by the non linear differential equation of Riccati type 
in the electrostatic analogy of exactly solvable pairing 
models.Ej Then we sketph the procedure originally em- 
ployed by Richardsona to obtain the thermodynamic 
limit of the BCS model. 

In the non-uniform case the electric field is 

E{z) := -fc(z) + _ + I ^ coth[p(e^ ~ z)\ (Bl) 

with k(z) as defined in the previous section. In order 
to have the field obeying a Riccati-type equation, it is 
crucial to map one of the isolated charges to infinity. 
This is done by the transformation 

gi := exp {2pei) ; Ca exp (2^6^). (B2) 

In the variables and Cq,, the Riccati equation is 



a = 1, . . . ,iV 



(A3) 



dz 



+ E\z) = f{z) + f{z) 



(B3) 
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where fo{z) := i/9j2j=i^j/[ljil3 - ^)] is a generating 
function of the eigenvalues of integrals of the motion of 
the model (P (see Eqs. (^) and 



m = 



QoiQo + 1) 



(B4) 



E 



kj{kj + 1) 



- 2- 



9j 



— - 2— 

2 



1/2. 



where Pq := 2(5o + ^-2N - 

For the isotropic limit p ^ the ECS rational case is 
recovered. The electrostatic field is here given by (see RS 
Eqs. (0)) 



N 



Eoiz) - ko{z) + --y2^— 



(B5) 



^ 1/2 ^ 



1 1/2 1 
= - + > — > 



In this case the Riccati equation reads 

+ i?2(z) = r(z) + /(z) 



where 



5 ~{ (^-^i 



(B6) 



(B7) 



where t(z) is the generating function of the eigenvalues of 
integrals of the motion of the BCS model Eqs. (p7|). The 
Richardson BA equations are the zeros of the solutions 
of (B^). The role of the Riccati equations in connection 
with the Gaudin models was also investigated in Ref.|3^. 

In both the rational or non-uniform cases it is evident 
that the Riccati equation plays an important role. This 
fact has been used by Richardson to derive the BCS 
theory from the expansion of the solution of the Ric- 
cati equation (in the rational case and kj = — l/2Vj) 
in powers of 1/N. In the following we briefly summa- 
rize the Richardson's physical arguments. The Richard- 
son equations, via the electrostatic analogy, can be seen 
as a self-consistent evaluation of the field of the fixed 
charges in the presence of the screening due to the mo- 
bile charges. The strategy followed by Richardson is to 
attempt to eliminate any reference to the mobile charges 
in the right hand side of the above equation. In fact, 
using the Eqs. ( [A^ ) the Riccati equation for the field can 
be recast in the following form: 



dEojz) 
dz 



+Ei{z) = y: 



1/2 



E 



1/2 



e/3 



(B8) 



In the last term on the right hand side there appears the 
field due to the mobile charges at the location of the fixed 
charges. The self-consistency condition, and the effective 
screening, enters here. Richardson notices that the field 
of the mobile charges may be written in the following way 



1 



ei - ep 2m e.^ - z 



(B9) 



with the curve C going around the singularities of E{z) 
due to the mobile charges. By using this, one gets a 
integro-differential equation of the Riccati type 



dz 



+ El{z) 



E 



1/2 



2^ 



:i=l 

Eoiz') 



1/2 

z - e., 



C fr^iz-£i)iz' -e,)' 



1 

f - 
9 

(BIO) 



The nice feature of the above equation is that any ref- 
erence to the mobile charges has disappeared. One may 
object that in the above equation there is still a reference 
to the mobile charges via the curve C entering the con- 
tour integral. In fact, by knowing the property that the 
field has to satisfy at infinity one may observe that the 
integral over a closed curve may be written as the sum 
of two curves: the first is C enclosing the singularities of 
the mobile charges, while the second enclosing the sin- 
gularities of the fixed charges. Furthermore to evaluate 
the contour integral along the curve extending at infinity, 
it is enough to know the first term in the multi-pole ex- 
pansion of the field. The final result is that the integral 
over C may be written in terms of the curve enclosing 
the singularities of the fixed charges and the behavior of 
the field at infinity only. The elimination of the mobile 
charges from the effective equation for the field is closely 
in the spirit of the Thomas-Fermi approach to the eval- 
uation of the effective potential due to the combined ef- 
fect of fixed charges (nucleus in atoms, ions in metals) 
and mobile charges (electrons). The integro-differential 
Riccati equation plays in this context the role of the Pois- 
son equation, which becomes the Thomas-Fermi equation 
once the charge of the mobile electron is written in terms 
of the potential. 

The next trick used by Richardson is to expand the 
complex electric field Eoiz) near the point at infinity, in 
powers of 1/z. This amounts to a multi-pole expansion. 
The constant electric field gives the zeroth order term. 
The other terms may be expanded upon using the geo- 
metric series. One easily gets 



E, 



E, 



(0) 


(1) 



1 



E, 



(2) _ 



1 

9 

N - -n 

2 

N 

Es'' 



(Bll) 



IE' 
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The monopole and dipole terms give then the number 
of pairs and the energy. To perform the thermodynamic 
hmit, we rescale the couphng constant as g ^ G/ft in 
order that the energy remains an extensive quantity in 
this hmit (see section III B). The field is expanded as 
Eq = En + El + Ei/pf + ... where the subscript indi- 
cated terms of order N, 1, and 1/A^, respectively. We 
are assuming that N and Q keep a fixed ratio. One may 
rewrite the Riccati equation order by order 



E% 



1Ej\jE\ — — 




(B12) 



where at each order in N 



27ri Jc ti- z 



(B13) 



Following the analogous lines presented in section III, 
Richardson wrote down the solution for E^q as 



E 



N 



1 



1 



(B14) 



with Z{z) = y/{ z ~ a)(z — b) having a branch cut be- 
tween the points a and b. To see how the BCS limit is 
recovered, one has to expand En in multiple's and insert 



the i?^^ back into (Bill). A direct calculation shows 
that 



E 



(0) 
N 



En ~ - 



E^ - - 



2 ^ Z(eO 
n 

l—l ^ ^ 

n 

2 ^ Z{e,) 



(B15) 



2 Q + ^ 1/ ,n2 



By setting a = A + iA and identifying the real and imag- 
inary parts of a with the chemical potential and the en- 
ergy gap, the first two equations of (B15) reproduce the 
gap equation and the condition for the chemical poten- 
tial, while the third gives the ground state energy. This 
completes the evaluation of the field in the leading order 
in 1/A^. To proceed further, one has to solve the second 
of (BIS ), which is complicated by the way the field Ei en- 
ters via the quantities Hu. Richardson showed, however, 
that there is an elegant way to circumvent the tackling 
of the integral equation. An important point to notice 
is that the singularities oi En exhaust all the charges, 
so that the correction Ei may have no poles at the po- 
sitio ns o f the fixed charges. Inspection of the second of 
Eq.( B12 ) shows that Ei has poles at the zeros of En. 
However, Ei cannot have poles because this would imply 



further charges. To avoid this one has to impose that the 
right hand side of the second equation of ( B12 ) has to 
vanish at the zeros of En- There are ft — 1 zeros of En 
and one gets — 1 linear equations for the Q unknown 
Hu. One more equation may be obtained by perform- 
ing a multi-pole expansion in the Riccati equation. The 
Riccati equation for l/z terms reads 



(0) p(l) 



2E^'"E, 







^2 ^2 



(B16) 



which may be rewritten as 



G 



By expanding in powers of 1 /N one gets 



(B17) 



N 



= 



(B18) 



The second of the above equation provides the fith lin- 
ear equation for the determination of the unknown Hu. 
These latter, o nce d etermined, may be inserted back into 
the second of ( B12 ) to obtain the correction Ei . 



APPENDIX C: ELECTRIC FIELD FROM 
CONTOUR INTEGRATION 

In this appendix we present the derivation of Eqs.([2^) 
and ( p5| ) exploiting standard results of complex analysis. 
Then we calculate the residue at infinity. Finally we dis- 
cuss some subtleties associated with the normalization 
condition for the charge density iy{z). 



1. The contour integral J dzE{z)/{z — z') 

Our first aim is to write the last term in Eq.(p2|) in 
terms of the ansatz field E defined by Eq.(p3}. To do 
this, we consider the contour G shown in Fig.(5|). From 
the generalized Cauchy theorem, we know that the sum 
of the contour integrations of E{z')/{z — z') inside G is 
equal to the negative residue of E{z') at infinity (i.e. the 
charge at infinity), which is zero in our case. 

The contours surrounding F can be contracted to an 
integration along the right-hand side minus an integra- 
tion along the left-hand side of F. Since the discontinuity 
of E(z) at z g F equals the charge density at z times 27ri , 
this gives 



1 

2Td 



dz 



E{z') 



Cr 



\dz 



a{z')-p{z') 



(CI) 
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that the leading order is 1/z^. 



, z J ^ — ' \ z 

j=0 ^ ' k=0 



where Cj are the coefficients of the expansion of 
nlo Vi^-an/z)il-bjz). 

Ij =Q+ + {-lyQ- + J dx x^f{x) . 

Equaling all powers z*^^ for i = Q,...,K to zero, we 
obtain the triangular set of linear homogeneous equations 
for /j : 

Cq/o = 

coh + Clio = 



cqIk 



crIo = 



FIG. 5: The integration contour L (full contour) and its defor- 
mations are shown. The density fj{z) of the "solution charges" 
is determined by evaluating once the inner part of L, and then 
the residue at infinity. The contour integration J E{z)/{z — z') 
enclosing the arcs of F extracts the charge density (j{z) on 
these arcs; f-^ \d z\a{z) / (z — z'). 



The integrals along the small half-circles at the end- 
points of r tend to zero. If we don't have an arc of F 
on L defined within the function S{z), we obtain 



=- 



27ri 



z — z' z — l/p 



S{-l/p)Q- 



+ I \dz 



cj{z')^p{z') 



S{x)ip{x) 



dx^O. 



(C2) 



This finally leads to 



aiz') 



\dz'\ 



Sil/p)Q+ S{-l/p)Q- 



z - 1/p 
S{x)(p{x) 
L\r z-x 



z + l/p 



dx 



p{x) 



dx 



Lnr 



(C3) 



By substituting Eq. (|C3|) in the Bethe equations Eq.(|22|), 
we find the conditions in Eqs.(p^. 

Note that a{z) and p{z) can have common support 
even though the screening function S{z) has no cut along 
L. Lnr here means the common support of (7{z) and 
p{z). 



2. Asymptotic behavior of E{z) at z ^ oo 
In order to find the remaining equations (|2^), we ex- 
pand the electric field in inverse powers of z, and impose 



which has the only solution (all Ci ^ 0) Iq = ■ ■ ■ = Ik 
0, which are the Eqs.(p5|). 



3. The normalization of a{z) from Eqs.{'^5\) 
We finally prove that the condition / \dz\ a{z) = 1/18 
contained within Eqs.(p5|). 

Contour integration of E(z) yields 

I \dz\(j{z) = [ dzE{z)+ / dzp{z) 

J 27ri J^np 

= ^S{l/p)Q+-S{-l/p)Q. 



dxS{x)ip{x)+ / dzp{z) 
L\r JLnr 



where we exploited the fact that Res {i?(z), oo} = 0, 
which is the equation in ( p5| ) for n = K if all the equa- 
tions for < n < X are fulfilled. 



APPENDIX D: GROUND STATE AND EXCITED 
STATES 

We argue that one should not include arcs with end- 
points on L. This is because the charges on such an 
ansatz arc have to be mobile. The charges on the real axis 
however find themselves arrested in the "cells" flanked by 
two adjacent fixed charges. Therefore the only possibil- 
ity for them to become mobile (in the thermodynamic 
limit) is to "escape" into the complex plane. This how- 
ever cao. happen only for solution charges in neighbored 
"cells"Ej. Zones of separated "cells" Lsep are character- 
ized by —p{x) := <j{x) — p{x) < 0; x G Lsep- We can 
contribute for them, just taking p as the density of fixed 
charges. The complementary region Lpairing '■= L\ Lsep 
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is thus characterized by a{x) — p{x) > at G = and we 
let p{x) := p(x). In the simplest case we have p(x) = p(x) 
everywhere and also the modulus of the total charge is 
p[x) in each point. This means a{x) — for x £ Lsep 
and distinct connected regions of neighbored "cells" with 
a{x) = 2p{x) for x e Lpairing at G = (see Fig.||). This 
most simple situation applies in particular to the ground 
state, which we have considered in the numerics. The 
general situation can be attacked replacing p{x) with p{x) 
and the filling v by 



d z|(t(z) 



G=0 



d xp{x) 



= V- i <lx[p{x) + p{x)] . 

That the above argument applies also to the general- 
ized class of BCS models discussed here, follows from 
the same structure of the equation after the transforma- 
tion. But it can equally be seen from the electrostatic 
analogy, since no two charges can penetrate each other 
due to the diverging forces when sitting upon each other, 
which must be zero instead due to the Richardson equa- 
tions. These infinite forces can only be overcome if ei- 
ther at least one of the external charges diverges (i.e. the 
case G=0) or if two solution charges together approach 
a fixed charge (in the presence of a degeneracy d, more 
than two solution charges are needed because the charge 
ratio \qmobiie/qfixed\ = 2/d is diminished by d). We come 
back to this discussion after having presented the solu- 



tion of the problem. Now we explain why the solution of 
Eq.(po[) corresponds to the ground state of the system, 
which implies that up to the crossing point of F with 
the Debye shell L all "cells" between fixed charges are 
occupied. This implies that there we have a — 2p and 
hence the total charge density is —p. This sign change 
of the charge density below the crossing point is already 
included in Eq.(pO), since we apparently replaced 



S{x) 



pix) 



(Dl) 



Notice that the complex square root changes sign at the 
crossing point Xc'- indeed we have 



Six) 



y^{x - A)2 + A2 X> Xc 



for 



-^{x - A)2 + A2 X <Xc 



This means that in Eq. (pO|) the charge density is —p{x) 
for X < Xc and p{x) for x > Xc, which correspond to the 
ground state. For general excited states, which might 
also include more than one arc, we must furnish the real 
square roots with the proper signs as indicated in Figj^. 
As long as all the arcs cross the real axis in the same in- 
terval of L. 



pairing 



where they are generated from, this pro- 
cedure can be implemented in the formalism described in 
this work. Otherwise the crossing points could be deter- 
mined numerically. 
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FIG. 6: Here it is indicated schematically how the signs of 
the real square roots have to be chosen in order to yield the 
sketched charge distribution. The function S{z) consists of 
complex square roots and has a sign change at each crossing 
point of r (dash-dotted grey curve) with the Debye shell (solid 
lines) . At the crossing point of the right-most part of F a sign 
change in the charge density {a = 2p) is intended and hence 
the sign change of 5* located there (dotted vertical line) is 
"shifted" to represent the desired sign change p — > —p. Grey 
parts of the Debye shell indicate a positive total charge density 
at G = and hence mobile charges; the black parts instead 
correspond to a negative total charge density at G = and 
hence separated "cells" . 
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In formulas ^ instead of the single particle energies the 
quantities rjj = Sj — gjj/2 + 2^^ Uij enter. We fixed gjj 
and Ujj such that rij = Sj 

We recall that in two dimensions the electric field generated 
by a point charge q having complex coordinate zq can be 
represented by the holomorphic function E — q/{z — zo)- 
The analogy is by no means unique. In fact, a general Moe- 
bius transformation p{z) — re^'^{z — loi)/{z — UJ2) uJi,ijJ2 G 
C preserves the algebraic structure of the Bethe equa- 
tions (|^, varying the locations of the charges. 
With ground state we mean that state connected with the 
ground state for G = 0. Up to some "critical" value of G 
this is the ground state of the system. 



